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Abstract. There is considerable uncertainty as to whether interannual variability in climate and 
terrestrial ecosystem production is sufficient to explain observed variation in atmospheric carbon 
content over the past 20-30 years. In this paper, we investigated the response of net C0 2 exchange 
in terrestrial ecosystems to interannual climate variability (1983 to 1988) using global satellite 
observations as drivers for the NASA-CASA (Carnegie- Ames-Stanford Approach) simulation 
model. This computer model of net ecosystem production (NEP) is calibrated for interannual 
simulations driven by monthly satellite vegetation index data (NDVI) from the NOAA Advanced 
Very High Resolution Radiometer (AVHRR) at 1 degree spatial resolution. Major results from 
NASA-CASA simulations suggest that from 1985 to 1988, the northern middle-latitude zone 
(between 30 and 60 degrees N) was the principal region driving progressive annual increases in 
global net primary production (NPP; i.e., the terrestrial biosphere sink for carbon). The average 
annual increase in NPP over this predominantly northern forest zone was on the order of +0.4 Pg 
(10 15 g) C per year. This increase resulted mainly from notable expansion of the growing season 
for plant carbon fixation toward the zonal latitude extremes, a pattern uniquely demonstrated in our 
regional visualization results. A net biosphere source flux of CO 2 in 1983-1984, coinciding with 
an El Nino event, was followed by a major recovery of global NEP in 1985 which lasted through 
1987 as a net carbon sink of between 0.4 and 2.6 Pg C per year. Analysis of model controls on 
NPP and soil heterotrophic CO 2 fluxes (Rh) suggests that regional warming in northern forests can 
enhance ecosystem production significantly. In seasonally dry tropical zones, periodic drought 
and temperature drying effects may carry over with at least a two-year lag time to adversely impact 
ecosystem production. These yearly patterns in our model-predicted NEP are consistent in 
magnitude with the estimated exchange of C0 2 by the terrestrial biosphere with the atmosphere, as 
determined by previous isotopic (5 13 C) deconvolution analysis. Ecosystem simulation results can 
help further target locations where net carbon sink fluxes have occurred in the past or may be 
verified in subsequent field studies. 
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Introduction 


Terrestrial biosphere response to interannual variability in climate may explain a significant 
portion of year-to-year variation in atmospheric carbon content observed over the past two decades 
(Keeling et al., 1989; Ciais et al., 1995; Keeling et al., 1996). It is plausible that terrestrial 
ecosystems introduce delay and damping mechanisms into the global carbon cycle that are 
detectable in the atmospheric CO 2 record. As evidence, fluctuations in atmospheric CO 2 associated 
with El Nino events appear to have lagged variations in surface air temperatures by one or two 
years (Keeling et al., 1995). However, there is a lack of direct evidence from field measurements 
to attribute these fluctuations with certainty to interannual variability in terrestrial ecosystem 
production. 

Satellite observations of vegetation greenness have been used to monitor the duration of the 
active growing season for terrestrial vegetation during the 1980s (Myneni et al., 1997). Longer 
growing seasons are apparent, particularly in areas of the northern high latitudes (between 45° N 
and 70° N), where notable warming has occurred in the spring. These satellite observations appear 
to be consistent with an increase in amplitude of the seasonal cycle of atmospheric CO 2 since the 

early 1970s. 


Nevertheless, greenness estimates from remote sensing data need to be combined with 
estimated controllers of ecosystem carbon fluxes in order to quantitatively assess how actual net 
CO 2 exchange has responded to fluctuations in climate (Running and Hunt, 1993; Goward and 
Prince, 1995). Ecosystem carbon models can take into account moisture stress and soil 
biogeochemical processes that may not be detected by satellite observations of vegetation 
greenness. The present study involves the global application of one such ecosystem carbon model 
that uses monthly 1° satellite daia to drive net terrestrial production estimates and process- level 
algorithms for soil CO 2 fluxes. 
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The NASA Ames ecosystem model version (Potter and Kiooster, 1997) of CASA 
(Carnegie- Ames-Stanford Approach) has been calibrated for non-equilibrium state simulations 
driven by satellite vegetation index data from the NOAA Advanced Very High Resolution 
Radiometer (AVHRR) over the years 1983-88, in an effort to capture some of the intrinsic 
ecosystem regulation of plant production and soil respiration by biophysical processes and 
interannual responses to anomalous events. 


Methods 

To make a direct assessment of terrestrial ecosystem exchange of carbon, we simulated 
global net primary production (NPP) and soil heterotrophic respiration (Rh) for the years 1983- 
1988 using the NASA-CASA Biosphere model version (Potter and Kiooster, 1997). As a unique 
feature of CASA, model predictions of NPP and R h have been validated previously against 
measured seasonal oscillations of atmospheric carbon dioxide (Potter et al., 1993, Denning, 1994), 
and against multi-year estimates of NPP from field stations and tree rings (Malmstrom et al., 

1997). 


The model is driven with a satellite greenness index from the NOAA Advanced Very High 
Resolution Radiometer (AVHRR). Calculation of monthly interannual NPP at all ice-free 
terrestrial grid cells is based on the concept of light-use efficiency. NPP is a product of cloud- 
corrected solar irradiance (S), fractional intercepted photosynthetically active radiation (FPAR) and 
a maximum light use efficiency term (e), modified by temperature (T) and moisture (W) stress 
scalars. 


NPP = S FPAR e T W 


( 1 ) 
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In the global or regional simulation mode, estimation of FPAR comes from a 
vegetation index derived from AVHRR FAS-NDVI data (Sellers et al., 1994), whereas the e term 
derives from calibration to previous field estimates of NPP. The T stress term is computed with 
reference to an AVHRR-derivation of optimal temperatures (Topt) for plant production (Potter et 
al., 1993). The T^t setting ranges from near 0° C in the Arctic to the middle thirties in low 

latitude deserts. 

Soil carbon cycling and Rj, flux components of the NASA-CASA model (Figure 1) are 
based on a compartmental structure comparable to the mechanistic CENTURY ecosystem model 
(Parton et al., 1992). First-order equations simulate loss of CO 2 from decomposing plant residue 
(metabolic and structural fractious of NPP) and microbes at the soil surface. Near-surface soil 
organic matter (SOM) pools are presumed to vary in carbon residence time and chemical 
composition. Active' (microbial biomass and labile substrates), slow (chemically protected), and 
passive (physically protected) fractions of the SOM are represented. 

Monthly 1° climate anomalies for air surface temperature and precipitation (Dai and Fung, 

1993) , together with long-term (30-year) mean values (Leemans and Cramer, 1990) and surface 
solar irradiance (Bishop and Rossow, 1991), were used to generate model drivers for the period 
1983-1988. Dai and Fung (1993) used, among other sources, the World Surface Station 
Climatology at NCAR. 

Model FPAR is derived using the Normalized Difference Vegetation Index (NDVI) from 
AVHRR (Sellers et al., 1994). The period 1983-1988 was used to eliminate the potential effects of 
satellite change-over and major volcanic eruptions. The monthly 1° FAS-NDVI product (Los et al. 

1994) includes Fourier smoothing algorithms (FA) and solar zenith (S) angle corrections to remove 
anomalous signals from the global NDVI data sets. Our calibrated global value (Potter et al., 1993) 
of e based on these climate and FAS-NDVI products is 0.57 g C MJ 1 PAR. 
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Complete AVHRR data sets for the 1980's have been produced from NOAA Global Area 
Coverage (GAC) Level IB data, and consist of reflectances and brightness temperatures derived 
from the five-channel cross-track: scanning AVHRR aboard the NOAA Polar Orbiter 'afternoon' 
satellites (NOAA-7, -9, and -1 1). Monthly composite NDVI data sets are designed to remove 
much of the contamination due to cloud cover present in the daily AVHRR data sets (Holben, 
1986). To generate a composite data set, eight to eleven consecutive days of data are combined, 
talcing the observation for each 5- to 8-km bin from the date with the highest NDVI value. Only 
data within 42 degree of nadir are used in the composite to minimize spatial distortion and bi- 
directional effect biases at the edge of a scan. A Rayleigh correction is calculated and applied using 
a standard radiative transfer equation and methodology, which follows the work of Gordon et al. 
(1988). 


As part of the Global Inventory Monitoring and Modeling Studies (GIMMS) program (Los 
et al., 1994) of NASA Goddard Space Flight Center, Sellers et al. (1994) describe Fourier 
smoothing algorithms (FA) and solar zenith (S) angle corrections for interannual AVHRR data sets 
to further remove anomalous NDVI signals from global 1° data sets (averaged from 8-km values) 
for the 1980's. The FA correction is intended to remove remaining cloud cover or aerosol 
interference. Settings for this FA correction include three temporal harmonics and a weighted 
Fourier transform, i.e. values which fall above the Fourier curve are given more weight than 
values below the curve. This assumes that higher NDVI values are more likely to be correct than 
low NDVI values which could occur during periods of cloud or smoke formation. FAS 
processing appears to eliminate artifactual problems present in other NDVI data sets. Specifically, 
FAS(IR) NDVI data sets have been shown to reproduce multi-year yield data and tree ring growth 
patterns from field studies world-wide (Malmstrom et al., 1997); the FAS(IR) NDVI data show 
minimal coiTelations with equatorial crossing times of the NOAA satellites, which suggests 
correction for orbital drifts and switches between satellites (e.g., NOAA 9 to NOAA 1 1). 
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FAS-NDVI data sets (scaled in units from 0 to 1000) for 1983-88 (Los et al., 1994) were 
made available for the present study. Unlike in earlier versions of the CASA-Biosphere model 
(Potter et al., 1993; Potter and Klooster, 1997), we avoided use in this NASA-CASA version of 
AVHRR-NDVI data sets with interpolation and reconstruction (IR) features added (Sellers et al., 
1994). The main justification for this decision is based on lack of verification of the NDVI 
reconstruction algorithm for the extent of tropical evergreen forest. 

To accommodate interannual variability in model drivers, we developed a protocol to 
initialize soil pools (water and carbon) with states that realistically represent near-term climate and 
vegetation greenness conditions preceding the 1983-88 time frame. The mean surface temperature 
over 1980-1982 (Dai and Fung, 1993) and the 30-year mean values for precipitation (Leemans and 
Cramer, 1990), along with mean monthly FAS-NDVI and solar irradiance from 1980s time series 
data sets, were found to most closely represent initialization drivers for the model to a near steady- 
state over a computational period of 300 years. We tested alternative driver combinations, but 
found that: (1) surface temperatures must reflect recent warming trends to produce near steady- state 
conditions for soil carbon (SLOW pool; Figure 1) typical for the late 1970s-early 1980s; (2) 
rainfall varies considerably from year-to-year, which justifies use of a long-term mean; and (3) 
NDVI data does not exist for the late 1970s-early 1980s, which justify use of a mean of existing 
data for the mid-1980s. 

The initial "near steady-state" pool sizes for water and modem carbon storage were used to 
start the interannual simulation run in the transient mode from 1983-1988 using actual (anomaly 
adjusted) monthly climate and FAS-NDVI drivers for this time period. Model state variables were 
allowed to change from month to month without the constraint of an equilibrium between NPP, 
litterfall transfer to the soil, heterotrophic soil CO 2 fluxes, or soil carbon turnover and storage. 
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Results and Discussion 

Our global model results identify several large regions as having experienced substantial 
interannual variability in NPP sink fluxes of carbon between 1984-88. In areas of Canada, 

Europe, and Russia, NPP increased from 1985-88 (Figure 2a). This is mainly as a result of 
warmer than long-term average spring temperatures and lower summer drought stress which 
promoted a northward expansion of plant production over the period (Figure 2b). Across the 
entire African Sahel, and in eastern and southern Africa and eastern Brazil, annual NPP showed a 
strong recovery from the severe drought effects of the 1983 El Nino event (Keeling et al., 1995). 

In areas of the north-central United States, North Africa, and Australia, mean estimated NPP 
decreased chiefly as a result of periodic shortages in precipitation. These model responses are the 
result of combined effects of NDVI and climate anomalies. Our regional analyses suggest that 
interannual fluctuations in the satellite greenness signal account for 30-80% of the variability in 
predicted NPP, with climate anomalies accounting for the remaining variability (also see analysis in 
Potter et al., submitted). This finding of NDVI accounting for more than 50% of the variation in 
modeled NPP is also consistent with the results of Malmstrom et al. (1997), who made interannual 
comparisons of CASA model predictions with agricultural yield data from sites around the world. 


The time series trend in predicted global CO 2 flux from the NASA-CASA model suggests 
that terrestrial ecosystem production was a net source for carbon in 1983 and 1984. Predicted 
NEP shifted to become a biosphere sink for carbon starting in 1985 at 2.6 Pg C yr l and 
continued, albeit at a diminishing rate, through 1988 (Figure 3). Higher plant production 
progressively adds organic material to regional soil pools which, with longer turnover times for 
carbon than in green vegetation, respond more gradually to warming in terms of CO 2 return to the 
atmosphere. The zone between 30° N and 60° N accounted for over 20% of the global net sink 
flux in 1987 and 100% of the carbon sink in 1988. These NASA-CASA results support other 
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accumulated evidence for a substantial terrestrial sink in the Northern Hemisphere (Keeling et al M 
1989; Ciais et al., 1995; Melillo et al., 1996 )., and help identify specific ecosystem-level 
processes for this inferred carbon storage. 

As a means of independent confirmation of NASA-CASA results, estimation of the 
terrestrial biosphere flux of carbon from isotopic (5 13 C) deconvolution analysis (Keeling et al., 
1995) appears most consistent with a net biosphere sink for CO 2 over the middle of the period 
(i.e., 1985-1987; Figure 3). A six to twelve month offset in peak CO 2 fluxes is evident in 
comparison of these two estimates. This offset might be explained solely by systematic biases 
cited in the deconvolution analysis (Keeling et al., 1995). We note furthermore that interannual 
variability in terrestrial biosphere fluxes of CO 2 from the 5 13 C deconvolution model method 
includes other (non-NEP) sources, such as biomass burning, which can be detected best with high 
resolution remote sensing (e.g., 30-m Landsat; Peddle et al., 1997). These biomass burning 
fluxes would help explain differences in the two model results such as those seen in 1983 and 
1987. Therefore, net carbon source fluxes are not strictly comparable in the two modeling 
approaches, a point demonstrated in global ecosystem analyses not driven by satellite data 
(Kindermann et al., 1996), but which also confirm the magnitude of NASA-CASA estimates. 


Four sites were selected to evaluate the model over a wide range of eco-climatic conditions 
for terrestrial carbon exchange (Figure 4). Net ecosystem production (NEP), defined as NPP-Rh, 
was used to characterize seasonal carbon exchange with the atmosphere. At all test sites, 
interannual fluctuation in NEP was mainly a function of variability in NPP fluxes rather than R h , 
which is consistent with results from previous global modeling studies of terrestrial ecosystem 
carbon fluxes (Kindermann et al., 1996). Predicted yearly R h shows substantially less year-to- 
year variation (coefficient of variation (CV) < 10%) over this time series than annual NPP (CV > 
10%), especially in the semi-arid grassland (West Africa) location, where the model predicts the 
highest interannual variability (CV > 25%) in NPP fluxes. At the Amazon forest location, 
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fluctuations in Rh track the seasonal patterns of monthly NPP with a lag time of just one to several 
months. This short lag results from the estimate of a relatively constant year-round supply of 
freshly deposited leaf material to soil decomposers in humid evergreen forests. Hence, the storage 
of relatively large pools of slowly decomposing litter and SOM carbon, which carry residence 
times of 2 to 25 years and generally account for 40-60% of annual Rh fluxes of CO 2 (Potter and 
Klooster, 1997) apparently buffers the interannual effect of climate variability on Rh and NEP over 

large areas. 

In contrast to the observed lagged effects in the Rh response, above ground fluctuations 
(10-50% from year-to-year) in mean annual FPAR have a direct and immediate effect on NPP sink 
fluxes of CO 2 . For example, at the northern boreal forest location in Canada, small but notable 
increases were detected in the 1° NDVI signal early in the growing seasons of 1984, 1985 and 
1988. This advanced spring greening is potentially associated with an early disappearance of snow 
cover, a pattern verified in comparatively high-resolution studies of the NDVI temporal signal over 
the years from 1981 to 1991 (Myneni et al., 1997). Warmer (than long-term average) air 
temperatures in the spring months can drive higher estimated plant transpiration and reduced 
surface water runoff. These changes in hydrology were combined with the observed greening 
trend in NDVI-FPAR for prediction of higher NPP by the NASA-CASA model. Several weeks of 
production were added during the early growing season. 

In addition, from 1985 to 1987, predicted drought stress during the summer decreased at 
the northern boreal forest site due to the combined effects of lower evaporative demand and 
sporadically higher rainfall. This can progressively extend the growing season in terms of NPP 
fluxes. The model result of 30-50% higher annual NPP over the years 1985-88 (relative to the 
"reference" year 1984, during which global climate anomalies were lowest) leads to the prediction 
of positive annual NEP at the middle- and high-latitude zones, typically in the range of +30 to +60 
g C m 2 yr 1 . These positive NEP values are the result of a small difference between two relatively 
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large numbers (NPP and Rh). However, each of these model predictions is qualitatively consistent 
with reported interannual controls on net ecosystem carbon exchange measured in tower studies of 
temperate forest sites (Goulden et al., 1996). 

Under warmer low-latitude (< 20° N or S) conditions, seasonal moisture deficits were seen 
to severely restrict monthly NPP at the western African location selected to represent semi-arid 
grasslands. The NDVI responded fairly strongly to interannual patterns of monthly rainfall, 
although predicted increases in peak seasonal greenness tended to lag positive precipitation 
anomalies by about a year (Figure 4; also see Potter et al., 1997). At the much wetter tropical 
forest location in the Amazon Basin (near Manaus), seasonal low NPP mainly showed the effects 
of notable rainfall deficits during the El Nino events of 1983 and 1987. Low rainfall in the 
Amazon reduces NDVI-FPAR estimates. Predicted soil Rh in the tropical moist forest location is 

less sensitive to these relatively severe drought events than NPP. 

In summary, we hypothesize that use of the satellite greenness index together with 
medium-term (decade long) climate fluctuations captures the most important responses of terrestrial 
ecosystems to interannual variability in biosphere-atmosphere exchange of C02- Analysis of 
NASA-CASA modeled map results (e.g., Figure 2a) can also locate new validation studies for net 
carbon sinks at the ecosystem scale. The spatial detail made possible through the use of satellite 
remote sensing data in ecosystem modeling uniquely facilitates this type of region-by-region, year- 
by-year assessment of potential C sinks. 

Although the theoretical basis for NDVI as a key land surface parameter related to plant 
carbon fixation has been spelled out previously in some detail (Running et al., 1996), much work 
remains with respect to combining optical remote sensing with radiative transfer modeling to 
enhance the accuracy of FPAR and leaf area index (LAI) estimates at scales from regional to global 
Assumptions and algorithms must be improved concerning leaf angle distribution, basic leaf optical 
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properties, and clumping of canopy components as a means to generate the improved relationships 
among satellite vegetation indices, LAI, and FPAR (Myneni et al., 1995). With improved 
algorithms in hand, a new generation of NASA Earth Observing System (EOS) platforms should 
provide decades of well-calibrated surface NDVI data for monitoring seasonal and interannual 
variability in terrestrial carbon sinks over regional to global scales (USGCRP, 1997). 

Results from this ecosystem modeling study, along with other studies that have used satellite 
observations to infer plant production, should be qualified nonetheless to point out that, while 
estimates may be indicative of interannual response of terrestrial carbon flux to variability in 
climate, their conclusions cannot be extrapolated directly to infer either relatively short-term (day- 
to-day) or long-term (several decades) response of the biosphere to global warming. The 
processes of ecosystem acclimation to higher surface temperatures, changes in daily surface 
hydrology, disturbances (e.g., fires), land use, and elevated atmospheric CO 2 levels must be 

included in subsequent modeling analyses that include potential geographic shifts occurring over 
several decades in plant functional types, physiological responses, and soil carbon turnover rates. 
The mechanisms controlling terrestrial NEP may interact differently under varying conditions of 
large-scale climate oscillation (i.e., changes in frequency and amplitude). The NDVI driver alone, 
while it has been proven useful as an indicator of potential past changes in plant greenness (Myneni 
et al., 1997), is not adequate to capture actual production changes of this extent and magnitude. 
Second-order soil contributions and feedbacks in lagged systems responses also require greater 
understanding and complementary modeling approaches. 
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Fig. 1. NASA-CASA soil sub-model. Litter and soil C and N transformations lead to substrates 
for trace gas production. Carbon pools are outlined in black and labeled with C-to-N ratios, C 
fluxes in solid arrows, CO 2 production in stippled arrows; Nitrogen pools in gray, N fluxes in 
gray arrows. Litter, microbe (MIC) and soil organic (SLOW and OLD) pools are shown. 

Structural (S) and metabolic (M) pools are shown for leaf and root litter. 

Fig. 2. Interannual results from the NASA-CASA model, a. Global map of mean terrestrial net 
primary production (NPP) over 1985-88, normalized by predicted NPP for the "reference" year 
1984. Yellow-red in the northern hemisphere temperate and high latitude zones indicate potential 
areas of carbon sink flux in response to warmer than long-term average spring temperatures and 
lower summer drought stress. Red in areas of the African Sahel and eastern Brazil indicate a 
recovery of annual NPP from the severe drought effects of the 1983 El Nino event. Areas in white 
indicate no detectable plant production or missing data. b. Predicted monthly pattern of NPP 
within the northern latitude zones over the model years 1985-1988, relative to the model 1984 NPP 
estimate. Higher peak seasonal production and the northward expansion of the growing season are 
evident over these years in comparison to 1984, during which annual NEP was close to zero for 
the region. Seasons are defined according to the months January, April, July, and October (J, A, 

J, O). 

Fig. 3. NASA-CASA model estimate (solid line) of global ecosystem carbon exchange with the 
atmosphere, compared to terrestrial biosphere flux of carbon recomputed from isotopic (8 13 C) 
deconvolution data (Keeling et al., 1995; dashed line). Running 12-month totals are plotted. 
Positive yearly mean values represent a net source flux from the biosphere to the atmosphere, 
whereas negative yearly values represent a net sink flux into the biosphere from the atmosphere. 

Fig. 4. Monthly model predictions for selected ecosystem sites. The top panels show estimated 
FPAR (dark solid line, unitless left vertical axis) and surface solar irradiance (dashed line, W nr 2 
on right vertical axis). The middle panels present estimated stress scalars for moisture (solid line, 
unitless) and temperature (dashed line, unitless); a value of one represents the lowest level of stress 
on NPP. The bottom panels show NPP (solid line) and Rh (dashed line), both as positive flux 
values in g C nr 2 . The initial model state is labeled as "INIT". Ecosystem type definitions 
(DeFries and Townshend, 1994) and geographic site locations are: northern boreal forest, 56° N 
126° W; temperate grassland, 53° N 9° E; semi-arid grassland, 16° N 12° W; tropical moist forest, 
2° S 59° W. 
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